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Abstract 

Background: Comparing biological time series data across different conditions, or different specimens, is a 
common but still challenging task. Algorithms aligning two time series represent a valuable tool for such 
comparisons. While many powerful computation tools for time series alignment have been developed, they do not 
provide significance estimates for time shift measurements. 

Results: Here, we present an extended version of the original DTW algorithm that allows us to determine the 
significance of time shift estimates in time series alignments, the DTW-Significance (DTW-S) algorithm. The DTW-S 
combines important properties of the original algorithm and other published time series alignment tools: DTW-S 
calculates the optimal alignment for each time point of each gene, it uses interpolated time points for time shift 
estimation, and it does not require alignment of the time-series end points. As a new feature, we implement a 
simulation procedure based on parameters estimated from real time series data, on a series-by-series basis, 
allowing us to determine the false positive rate (FPR) and the significance of the estimated time shift values. We 
assess the performance of our method using simulation data and real expression time series from two published 
primate brain expression datasets. Our results show that this method can provide accurate and robust time shift 
estimates for each time point on a gene-by-gene basis. Using these estimates, we are able to uncover novel 
features of the biological processes underlying human brain development and maturation. 

Conclusions: The DTW-S provides a convenient tool for calculating accurate and robust time shift estimates at 
each time point for each gene, based on time series data. The estimates can be used to uncover novel biological 
features of the system being studied. The DTW-S is freely available as an R package TimeShift at http://www.picb.ac. 
cn/Comparative/data.html. 



Background 

Comparing and characterizing temporal changes in gene 
expression is a routine method to elucidate functional 
features of biological processes [1,2]. If gene expression 
levels are measured across time in different individuals, 
or under different conditions, it may be desirable to 
compare the timing of gene expression changes between 
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two or more time series. As the number of large-scale 
time series experiments increases, so does the need for 
detailed and comprehensive comparisons between them. 
Such analyses can reveal dynamic features of the system, 
such as acceleration or delay in temporal expression 
changes caused by conditional perturbations, or indivi- 
dual or evolutionary differences. One example of such 
an analysis was the study of developmental gene expres- 
sion patterns in two insect species: mosquito {Anopheles 
gambiae) and fruit fly (Drosophila melanogaster) [3]. 
This analysis, based on the temporal alignment of the 
gene expression profiles between the two species, 
revealed distinct groups of genes involved in functions 
known to be specific to one of the species, including 
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delay in mosquito cuticle synthesis, as well as faster 
maternal-zygotic transition in flies. Other recent applica- 
tions of temporal alignment algorithms include: reten- 
tion time alignment of complex LC-MS data in 
proteomics and metabolomics studies [4,5], and align- 
ment of transcriptome changes in two types of the car- 
diac hypertrophy [6]. 

The methods employed to address such temporal 
alignment problems are commonly based on the 
dynamic time warping (DTW) alignment method. DTW 
was originally developed for dealing with speech recog- 
nition problems [7,8]. It utiUses dynamic programming 
to find an optimal alignment with respect to a given 
scoring function. The application of DTW to gene 
expression data was pioneered by Aach and Church [9] 
and has been further developed by other groups [10-14]. 

In this study, we tested whether the temporal align- 
ment of ontogenetic expression profiles can reveal novel 
biological features of human brain development. To do 
so, we used published gene expression time series from 
the human, chimpanzee and rhesus macaque prefrontal 
cortex [15], as well as human and rhesus macaque cere- 
bellum [15]. Human development is known to proceed 
more slowly than development of chimpanzees and 
macaques. This trend is reflected in longer gestation 
time (humans - 280 days, chimpanzees - 220 days, rhe- 
sus macaques - 165 days), later age of sexual maturity 
(in humans - 13 years and in the wild; in chimpanzees - 
8 years and in rhesus macaques - 6 years), and longer 
maximal lifespan (humans - 100 years, chimpanzees - 60 
years, rhesus macaques - 40 years) [16,17]. Further, tran- 
scriptome changes during prefrontal cortex development 
tend to be delayed in humans compared to chimpanzees 
[15]. On the phenotypic level, delay in human develop- 
ment with respect to the closely related primate species, 
such as chimpanzees, was long hypothesized to play an 
important role in evolution of human cognition [18,19]. 
Thus, investigation of the temporal differences between 
developmental gene expression profiles in the human, 
chimpanzee and rhesus macaque brains using DTW- 
based alignment method might shed new light on the 
biological processes underlying human cognitive 
development. 

The common features of our algorithm and other 
alignment methods [11] include: (a) use of interpolated 
expression values for two aligned time series based on 
fitted higher order spline curves, and (b) use of one of 
the time series as a template for aligning the second 
one, using a resampling procedure based on Euclidian 
distance minimization. Our alignment approach, in 
common with other methods, cannot be applied to 
genes that have substantially different temporal expres- 
sion behaviour in different time series. In such cases, 
the difference between two time series is likely to 



involve factors other than a change in expression tim- 
ing and, therefore, time shift estimation may be 
inappropriate. 

One of the advantages of our method, compared to 
the original DTW algorithm, is that it follows the 
approach first suggested by Smith et al, [11] in lifting 
the requirement to align the ends of the two time series. 
This is a small but important modification as for many 
naturally occurring processes the timing of a process' 
commencement and termination, as well as the timing 
of corresponding points across different time series, are 
not known. This improvement to DTW, introduced in 
our method, is not based on local alignment as in Smith 
et al, [11], nor based on segment alignment as in Smith 
et al, [14], but on a dynamic programming-based align- 
ment procedure adopted during global alignment. 

The novel aspects of our algorithm include it's ability 
to not only estimate the optimal alignment for each 
time point of each gene, but also determine the reliabil- 
ity of this estimate. Specifically, we implement a multi- 
ple simulation procedure based on real data parameters 
estimated on a gene-by-gene basis. This allows us to 
assess the significance of estimated time shift values for 
each time point of each gene. Furthermore, we are able 
to determine the global significance level for each gene 
based on false positive rate (FPR) estimates. The signifi- 
cance estimation is based on implementation of our 
dynamic programming-based algorithm, which is cap- 
able of both simulating thousands of expression profiles, 
and determining simulated time shift values, at fast 
speed. This is a unique feature of our algorithm com- 
pared to other similar alignment methods. 

Our method allows the user to select genes with time 
shift estimates determined to the desired degree of con- 
fidence. The time shift profiles that are obtained, can be 
further used to group genes according to their align- 
ment properties, rather than their original expression 
profiles. When we applied our method to published 
gene expression time series of human and chimpanzee 
prefrontal cortex development, we were able to deter- 
mine additional biological information by separating 
genes showing human-chimpanzee expression diver- 
gence into distinct groups, based on their time shift 
properties. Furthermore, when we applied our method 
to published gene expression time series of human and 
rhesus macaque prefrontal cortex and cerebellum devel- 
opment, we found that genes expressed in the gray and 
white matter of the brain show synchronized time shifts 
between the two species, despite differences in expres- 
sion profiles between the brain regions. These examples 
demonstrate the utility of the DTW-S algorithm, which, 
by implementing the significance level of the gene-by- 
gene variable time shift estimates, allows us to obtain 
novel biological insights from time series data. 
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Methods 

The DTW-Significance algorithm (DTW-S) 

There were three main motivations to advance the exist- 
ing DTW algorithm: (i) to find the best alignment 
between two time series, (ii) to estimate the value of 
time shift between these series at each time point and, 
(iii) to estimate probability to observe a given time shift 
at each time point by chance, when actual time shift 
between two time series is zero. 

DTW-S works as follows: we consider two series, x 
and J, with the number of time points m and n respec- 
tively: {xi,,„yX^) and (yi,...,Jn)> where n and m can be dif- 
ferent. Since expression time series tend to be noisy 
(due to both technical and biological sources) and sam- 
ples' ages may not be uniform across the age-range 
(human/primate samples are often collected opportunis- 
tically), we use cubic spline interpolation, as implemen- 
ted in the smooth, spline R function, to construct the 
age-expression trajectory for each gene within a sample 
group. This is similar to the approach of Smith et aL 
[11]. To avoid overfitting, a generalized cross-validation 
procedure using a suitable smoothing parameter is 
implemented in the smooth, spline R function [20]. The 
degrees of freedom parameter for spline curve fitting of 
the actual biological data was set to four. This para- 
meter provided a good data fit for the actual number of 
data measurements (39 for humans and 14 for chimpan- 
zees), determined by visual inspection of individual gene 
expression profiles. We subsequently reconstruct unob- 
served expression levels at uniformly distributed time 
points. This, in turn, enables us to estimate the time 
shift between two time series with greater precision. We 
denote the interpolated time series by X and Y: (X^,..., 
Xm) and (n,..., Y^), 

Our dynamic programming-based algorithm aims to 
find an optimal alignment between two time series (X 
and Y with number of time points M and N respec- 
tively, M>N suppose we align Y to X), while allowing 
two conditions: 

(1) Multiple mapping, i,e, one time point in one time 
series may correspond to many time points in the other 
time series, and vice versa. 

(2) Balanced alignment. We require the same number 
of time points from the two time series, to be selected 
for alignment. That is, we will have exactly (M-N) 
omitted points, assuming M>N, 

With the above two constraints, the goal is to find a 
pathway with least weight from the bottom row to the 
top row, with the nodes in the path covering exactly N 
columns. This additional information that covers exactly 
N columns cannot be kept in the subproblem since we 
will query the subproblem with less columns. Therefore, 
we define a warping matrix as Mat[A^] [M] [M-A^], and 
the element of Mat as Mat[i] [j] [k]. This is the best 



solution to align the first i points in Y to the first j 
points in X, and only (j-k) points in X will be aligned to 
the first i points in Y. 

In the above procedure, for each cell, time complex- 
ity of 0{M) is required. By introducing an indirect 
recursion and maintaining two matrices, we can 
decrease the time complexity from 0{M) to 0(1) for 
each cell, and express the time complexity for the 
algorithm as 0(M x N x (M-N)), Assuming N, M are 
the same order, the time complexity for our algorithm 
is 0(M^3). In comparison, if 0{M) is required for each 
cell, the time complexity will be 0(M^4). For example, 
running our DTW-S algorithm on 1,000 genes with M 
= 40 and N = 20, on a machine with a single core pro- 
cessor (1.83 GHz speed and 2 GB memory), takes only 
12 seconds. 

We use an example gene expression profile (VTAl) of 
human and chimpanzee brain maturation to illustrate 
our procedure [15]. The dataset contains gene expres- 
sion time series from 39 humans (0-47 years of age) and 
14 chimpanzees (0-44 years of age). The ages of sampled 
individuals are distributed non-uniformly between birth 
and maturity (Additional file 1, Figure SIA). First, we 
estimate gene expression trajectories using splines 
(cubic polynomial interpolations with four degrees of 
freedom). Next, we interpolate gene expression values at 
40 time points within human time series (X) and 20 
time points within chimpanzee time series (Y), distribu- 
ted uniformly over the sampled time intervals {M = 40, 
N = 20). An example of this procedure is illustrated in 
Figure lA. 

Next, we align interpolated values of the Y time series to 
the X time series using the dynamic programming-based 
algorithm as described above. As a result, among all 
possible alignment combinations that preserve the order 
of the time series, the one with the minimal Euclidean 
distance between N time points subset from a total of 
M time points of X time series, and N time points from 
Y time series, is selected as the optimal alignment 
between X and Y. An alignment example is shown in 
Figure IB. 

Based on this optimal alignment, we further calculate 
the time shift between the two time series as the differ- 
ence between time of point within Y and the subset of 
X. Specifically, if Yi is aligned to Xj, shift(Yi) = time(Xj)- 
time(Yi). If Y^ is aligned to (Xj,..., Xj+q), shift(Yi) = mean 
(time(Xj),..., time(Xj+q))- time(Yi). We illustrate the opti- 
mized shift for an example gene in Figure IC. 

Since our algorithm does not suit the analysis of genes 
that have substantially different temporal expression 
behavior in different time series (as such genes might 
involve factors other than change in expression timing), 
we recommend that users filter the tested genes by a 
correlation test. 
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Figure 1 Illustration of time series alignment using DTW-S. (A) The x-axis shows human (dark yellow) and chimpanzee (yellow) age on a 
log2-scale, including conception time following [15]. The y-axis shows expression of an example gene {VTAl) on a log2-scale. The triangles and 
the diamonds represent real expression values. The curves represent the spline curves. (B) The alignment results for the DTW-S procedure. The x- 
axis represents the rank of ages for interpolated points. The circles on the curves represent the interpolated expression values. The blue lines 
show alignment matches between the two interpolated time series calculated by the DW-S algorithm. (C) The optimized time shift estimates 
for chimpanzee ages based on the expression pattern of VTAl, after aligning the chimpanzee time series to the human time series. (D) The y- 
axis shows the expression level after shifting chimpanzee expression time series by the magnitude of the time shifts estimated by the DW-S 
algorithm (shown in C). 



The warping algorithm is written in C and R lan- 
guages and is available within the R package TimeShift 
at http://www.picb.ac.cn/Comparative/data.html. 

Estimating the significance of time shifts and the false 
positive rate 

A novel feature of our algorithm is that it calculates the 
significance of the time shift estimate per gene (with the 
null hypothesis of zero time shift) and determines the 
FPR. For this, we designed the following simulation pro- 
cedure: suppose that we have Ng genes, each of which 
contains two gene expression series x and y, with length 
m and n respectively. We apply the above-described 
DTW-S steps to these A^^ genes and calculate the time 
shift gene-by-gene for each time point. We call these 
estimates "real shifts". Then, for each gene, we construct 
a background zero-time-shift distribution. Specifically, 
we simulate two time series, based on one of the two 
original time series, with zero time shift between them 
and estimate "simulated shift" values 1,000 times. Next, 
we create 50 more simulated time series with zero time 
shift, and construct 1,000 background zero-time-shift 
simulations for each of them. Finally we compare differ- 
ences between the real shift and the background zero- 
time-shift distribution, and 50 simulated zero-time-shift 
time series, to their background zero-time-shift distribu- 
tions. We then estimate the final significance of the 
time shift, based on the difference between these two 
differences. 

In detail, the procedure runs as follows: 

1. For each gene, build a spline model called Smodel 
to fit one of the time series (either x or y) from the 



original dataset and calculate the residuals. Next, con- 
struct a model error distribution (Nerror) as a normal 
distribution, with the mean and the variance equal to 
the mean and the variance of the residuals. 

2. Simulate two time series with zero shift by model- 
ing the expression levels for the real time points of the 
X and y time series, using the same Smodel for both and 
adding errors based on the Nerror distribution. Repeat 
this step 1,000 times. 

3. Apply DTW-S to the simulated time series to 
obtain shift estimates for each time point of each gene. 
For each time point, count the proportion of simula- 
tions with the absolute value of shift greater than or 
equal to the absolute value of real shift. Calculate a p 
-value for randomly observing a shift as large as that 
observed. Then, count the number of significant time 
shift points per gene (in the analysis we used the nom- 
inal cut-off p < 0.05 for choosing significant time shift 
points). We term the distribution of the number of sig- 
nificant time points across all tested genes, the "real" 
distribution. 

4. For each gene, generate 50 extra pairs of simulated 
time series with zero shift following steps (1) and (2), 
and estimate the time shift by DTW-S. For each of 
these 50 zero-time-shift time series, compare the time 
shifts to 1,000 simulated shifts in the same way as com- 
paring the real shift to 1,000 simulated shifts in step (3). 
So, for each of the 50 simulations, and for each time 
point of each gene, we can determine how frequently 
zero-time-shift time series yielded significant time shift 
estimates compared to the 1,000 simulated shift distri- 
butions (false positives). Note that these false positive 
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estimates are based on time series simulated using real 
data parameters on a gene-by-gene basis. Finally, count 
the number of significant false positive time shift points 
per gene, for each of the 50 zero-time-shift time series, 
and build a corresponding distribution of the number of 
significant false positive points across genes, called a 
"null" distribution. 

5. For each of the significant time shift points per 
gene, we calculate the FPR as the proportion of genes 
with an equal or greater number of significant time shift 
points per gene in the "null" distribution, compared to 
the "real" distribution. Based on this comparison, the 
minimal number of significant time shift points per 
gene can be estimated at the desired FPR level. A gene 
is classified as having a significant shift if it contains the 
number of significant time shift points equal to or 
greater than this minimal number. 

Preparation of simulated datasets with known time shift 

1. We generated a time series x with length m from a 
given type of model (x = f(t)). We used three types of 
function: sine, linear and polynomial (described below). 
To each of the n time points, we added random error 
sampled from a normal distribution. This gave us a first 
expression time series with time t 

2. We generated a second expression time series y 
with the same length n, using the same model (y =f(t)), 
and the same error distribution, as the first time series. 
However, for this series we shifted the time t in the 
model calculation by a shift of known magnitude: time 
shift = Aj,..., A„, where A denotes shifts at each time 
point, which can themselves be a function of t. 

3. We applied DTW-S on x and y, by first interpolat- 
ing M and N uniformly distributed time points across x 
and y, to generate X and Y respectively (with M>N). 
Then we estimated the time shift of aligning Y to X by 
applying DTW-S. 

4. The entire procedure was repeated 1,000 times, 
starting from the simulation of x and y. We finally com- 
pared the obtained time shift estimates to the original 
set values. 

Results and Discussion 

Testing DTW-S using heterochrony simulations 

To test the efficiency of time shift estimation by DTW- 
S, we first applied the algorithm to a number of simu- 
lated datasets with known heterochrony. 
Sine function 

We tested whether our method could be applied to peri- 
odic expression changes, such as expression changes 
associated with the cell cycle or circadian rhythm. 
Across the trials, we used the same error distribution 
parameters and time shift, while we varied the three 
parameters n, N, and M (these parameters reflect 



numbers of actual and interpolated time points, see 
Methods for Details). For simplicity, we assumed num- 
bers of actual and interpolated time points to be equal 
{n = N), The expression changes with time were mod- 
eled as 3/ = sin(TT''25/t)+s, with error e~A/'(0,0.3) and fixed 
time shift A = 5 for all time points. Using these error 
distribution parameters, residual variance constituted 
approximately 15% of the total variance in the simulated 
datasets. This percentage of error-related variance is 
within the range of technical and biological variance, 
relative to factors such as age and species identity, 
observed in actual microarray experiments {e.g, [15]). 
These results demonstrate that DTW-S estimates are 
robust with respect to the numbers of original data 
points and interpolated time points. Furthermore, 
DTW-S could estimate the true time shift value (Figure 
2). As expected, variation among time shift estimates 
decreased with increasing number of time points. 
Linear and polynomial functions 

We tested DTW-S on non-periodic gene expression 
changes, such as changes described by linear and poly- 
nomial curves. Most organismal processes, e.g, develop- 
ment, aging, or response to stress, involve such non- 
periodic patterns. Specifically, we tested the algorithm's 
ability to determine constant or variable time shifts, as 
well as it's performance at different levels of error var- 
iance (s) and with a different sample size. We tested the 
algorithm using the following parameters: n = 20, N = 
20 and M = 40. 

First, we applied DTW-S on simulated time series 
based on linear and polynomial curves with a constant 
shift. We used the function y = a+bt+dt^+z, with a = 
(0,1), b = (-1,0,1), d = (-1,0,1), to simulate expression 
time series. In addition, we used a fixed time shift A = 2 
for all time points. We set the error distribution para- 
meters at -10% or -20% of the total variance, Le, 10% 
or 20% of the total variance is attributable to error. This 
proportion is within the error variance range observed 
in actual microarray experiments. Under both error 
rates and for both linear and polynomial curves, our 
method yielded accurate and robust time shift estimates 
(Additional file 1, Figures S2, S3). 

Second, we tested the ability of our method to deter- 
mine variable time shifts between two time series. In a 
biological system, it can be reasonably assumed that the 
difference in timing of expression changes between two 
processes might increase or decrease gradually. We 
therefore used four sets of variant time shift (the linear 
shifts: CI and C2, and the polynomial shifts: C3 and C4) 
to simulate expression time series, two linear and two 
quadratic (y = 3+2^-^^+s and y = 3-2^+^^+s), with error 
s~A/^(0,4), such that residual variance constituted -20% 
of the total variance. Our results showed that under 
these conditions DTW-S can effectively identify the 
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Figure 2 Estimating constant time shift between two simulated time series. (A) The data shown are expression levels of time series for 
species 1 (spl, purple) and species 2 (sp2, light blue), modelled using the sine function with different sets of N and M parameters. We use the 
DVM-S procedure to align the sp2 time series to spl, and calculate time shift values for the ages of sp2. The parameters N and M shown above 
the panels represent the numbers of interpolated time points for sp2 and spl time series, respectively. The number of the original time points n 
= N. (B) Time shift estimates calculated by the DTW-S procedure. Time shift estimates are based on 1,000 simulations of the expression time 
series shown in (A). The gray lines show the actual time shift values used in for data modelling. The parameter r represents the proportion of 
the total variance introduced by simulated error (s). 



modeled variable time shift between two time series 
(Figure 3). Specifically, in simulations for non-linear 
shift patterns (C3 and C4 in Figure 3), the predicted 
time shifts were significantly better described by the ori- 
ginal polynomial shift models, than they were by linear 
models (ANOVA p < 0.01). 

Next, we assigned the same four sets of variant time 
shift (CI, C2, C3, and C4) to time series, modeled using 
combinations of the expression profile parameters a, b, 
d {Le, linear and polynomial curves), with -20% residual 
variance. We found that, irrespective of the expression 
curve type and the time shift type, our method robustly 
identified the original shift profile between time series 
(Additional file 1, Figure S4). 

Lastly, we tested our method's performance under 
30% residual variance. Overall, DTW-S performed well 
at this relatively high error rate. Compared to the origi- 
nal time shift, the estimated time shifts showed higher 
variance at the middle of the age distribution, and biases 
towards positive and negative time shift estimates at the 
beginning and end of the age distribution respectively 
(Additional file 1, Figure S5). We note that, under even 
higher levels of noise, DTW-S might not be able to dis- 
tinguish a non-linear time shift trajectory from a linear 
one. In these cases, it might still be possible to 



accurately estimate the magnitude and type of shift by 
averaging across multiple genes to minimize noise. 

Application to primate brain development 
Human and chimpanzee cortical development 

We first applied the DTW-S algorithm to investigate the 
heterochrony patterns between the brain developmental 
profiles of humans and chimpanzees. Prevalence of a 
specific type of heterochrony, delayed development in 
humans (or neoteny), was previously reported in human 
and chimpanzee brains [15]. This result was based on 
expression time series measured in the dorsolateral pre- 
frontal cortex in 39 humans (0-47 years), 14 chimpan- 
zees (0-44 years) and one outgroup species: 9 rhesus 
macaques (1-18 years). For all three species, the dataset 
covered most of the postnatal ontogenesis and matura- 
tion, from birth until adulthood (Additional file 1, Fig- 
ure SIA). In this study, we re-examined this work using 
DTW-S. Our hypothesis was that expression hetero- 
chrony patterns might vary among genes, even among 
those that show the same type of heterochrony, with 
regard to the timing of the shift, as well as shift magni- 
tude. Following the original analysis, we first selected 
1183 genes satisfying the following criteria: (1) signifi- 
cant expression changes with age (using polynomial 
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Figure 3 Estimating simulated time shifts with DTW-S. The leftmost columns show the distributions of the four sets of variable time shift 
(CI, C2, C3 and C4). In panels CI and C2, the simulated time shift is increasing or decreasing linearly, while in panels C3 and C4, the time shift 
follows a quadratic trajectory. The second column from the left shows the corresponding two expression time series, spl (purple) and sp2 (light 
blue). These are simulated based on the function y = 1+t+t^+s (y is the expression level, t is time and s shows error), as well as on the time 
shifts following the patterns CI to C4. The third column from the left shows the two expression time series after shifting sp2 to spl by the 
magnitude of the shifts estimated by the DTW-S. The rightmost columns display the corresponding time shift values estimated by the OVM-S 
algorithm in 1,000 simulations of spl and sp2 time series, with random error s~A/(0, 4). 



regression models, at F-test p < 0.05), (2) significant 
expression difference between humans and chimpanzees 
(using analysis of covariance (ANCOVA), at F-test p < 
0.05) and, (3) significant positive correlation between 
human and chimpanzee expression profiles (Pearson 
correlation p < 0.05 and r > 0). We then applied the 



DTW-S algorithm to these 1183 genes to identify 
expression trajectories with significant time shift 
between the two species. Here, we aligned the chimpan- 
zee time series to the human age-scale. For each gene, 
we set the p-yslue cutoff for the significant time shift 
points at 5%, and the gene significance cutoff as more 
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than one third of the gene's time points having signifi- 
cant time shift (Methods). At these cutoffs, 482 out of 
the 1183 genes tested had significant time shift, and the 
FPR was estimated at 10.7% (Figure 4A). 



We assigned these genes to one of four phylo-ontoge- 
netic categories, using our previously reported classifica- 
tions [15]: (a) human-specific neotenic expression (here 
neoteny denotes a slow-down or delay in development). 
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Figure 4 Human-chimpanzee expression heterochrony analyzed by DTW-S. (A) Numbers of significant time shift points per gene for 1183 
genes tested for heterochrony in the primate brain expression dataset. The x-axis shows the number of interpolated time points (out of a total 
of 20) with significant time shift (p < 0.05) estimated using DTW-S (blue), as well as the mean time points in 50 simulated datasets (light blue). 
The simulation is based on the human expression time series (see Methods). In this case, the FPR is -10% for genes with >7 points showing 
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significant time sliift (denoted by tine liglit pinl< sliaded area). In tiiis analysis, we have aligned chimpanzee expression time series to human time 
series. (B) Proportions of genes assigned to four phylo-ontogenetic categories. The labels are the four phylo-ontogenetic categories (see Results). 
The colors indicate gene classification into the phylo-ontogenetic categories stated in the original publication [15], carried out using the DW-S 
algorithm by aligning chimpanzee expression time series to human time series, or aligning human expression time series to chimpanzee time 
series. (C) Patterns of time shift between human and chimpanzee time series for genes with human-specific neotenic expression, calculated by 
aligning chimpanzee expression time series to human time series. The numbers above each graph represent numbers of genes within the 
cluster. The upper row panels show time shift estimates for mean shift based on all genes in a cluster (red lines) and the distribution of the 
shifts (arrows) for individual genes within clusters (5%-95%). Positive values of time shift indicate that human ages are mapped to younger ages 
in the chimpanzee time series, indicating slower/delayed human development. Time shift shown on the y-axis is calculated as the base-two 
logarithm ratio between chimpanzee age and human age at aligned expression time points. The bottom row panels show the relationship 
between human and chimpanzee ages, based on the calculated time shift estimates shown on the upper row (yellow lines). Shaded areas show 
the human age range corresponding to the maximum shift interval. (D) Panels show the relationship between human and chimpanzee ages, 
based on all the calculated significant time shifts between human and chimpanzee time series for genes with human-specific neotenic 
expression. The genes are identified as neotenic either by aligning chimpanzee to human time series, or by aligning human to chimpanzee time 
series. The numbers above each graph represent numbers of genes within the cluster. (E) Properties of genes in three clusters defined based on 
their human-chimpanzee time shift patterns as shown in (C). The left panel shows the proportion of genes in each of the three clusters that 
overlap with genes that are expressed preferentially in a particular cell type: neurons, oligodendrocyte precursor cells (OPC), oligodendrocytes, 
astrocytes or non-specific. The right panel shows the proportions of genes in each of the three clusters that overlap with genes that are 
expressed preferentially in gray matter, white matter or non-specific. 

V J 



(b) human-specific accelerated expression, (c) chimpan- 
zee-specific neotenic expression and, (d) chimpanzee- 
specific accelerated expression (Additional file 1, Figure 
S6). To do so, we first determined the direction of the 
time shift (acceleration, neoteny, or none) between the 
human and chimpanzee time series, for all 482 genes 
with significant time shift. If at least 70% of time points 
showed consistent shift direction, we classified the 
related gene as accelerated or neotenic. We then deter- 
mined whether expression change has happened to the 
human or to the chimpanzee evolutionary lineage, using 
the rhesus macaque time series as the outgroup (Meth- 
ods). For all 482 genes, we could assign 337 (at lineage 
assignment test p < 0.05), or 267 (at lineage assignment 
test p < 0.01), to one of the four phylo-ontogenetic cate- 
gories. This analysis yielded a clear excess of genes in 
the human-specific neotenic expression category (Figure 
4B), which adds support to our previously published 
results [15]. 

To test the robustness of our method, we reversed 
the entire time shift identification procedure by align- 
ing the human time series to the chimpanzee trajec- 
tory. This yielded 642 genes with significant shift (FPR 
= 8.7%, /7-value < 0.05, number of significant time 
shift points per gene = 7). The difference in the num- 
bers of genes with significant time shift compared to 
the previous procedure, is most likely related to sam- 
ple size differences between the chimpanzee and the 
human time series (Methods). Despite this variation, 
the two sets of genes identified in the two procedures 
largely overlap: 369 between 482 and 642 genes 
(Fisher exact test: p < 0.0001). The proportions of 
genes assigned to the four categories were also highly 
consistent (Figure 4B). 



In the previous study, any differences in time shift 
magnitude among the heterochronic genes, were 
ignored. Furthermore, the time shift was assumed to 
remain constant across all time points [15]. To test 
these assumptions, we analyzed the time shift patterns 
estimated by DTW-S for 118 human-specific neotenic 
genes. Using hierarchical clustering we found that, 
instead of forming a single pattern, the time shift pro- 
files of these genes fell into three distinct time shift clus- 
ters (Additional file 1, Figure S7). Clustering time shift 
profiles using the k-means algorithm also produced con- 
sistent segregation into three clusters. In 828 out of 
1,000 iterations, the k-means algorithm yielded the same 
distribution of time shift patterns across the three clus- 
ters. The three clusters also showed marked differences 
in both the timing and the ampUtude of the time shift 
(Figure 4C). Specifically, all three clusters showed time 
shift peaks at 20-30 years of age in humans, but at dif- 
ferent ages in the chimpanzees: 1 year (Cluster2), 2-4 
years (Cluster3), and 4-6 years (Clusterl) (Figure 4C). 

To test the robustness of the time shift estimates for 
the three clusters, we combined the estimates obtained 
in the chimpanzee-to-human alignment, and in the 
human-to-chimpanzee alignment, and repeated the clus- 
tering using the k-means algorithm. This approach also 
led to robust segregation into three clusters. In addition, 
for the 170 genes with significant time shift estimates in 
at least one of the two alignments, both direct and 
reverse alignments produced consistent time shift esti- 
mates between human and chimpanzee expression pro- 
files within each of the three clusters (Figure 4D). 

Finally, we investigated whether the various patterns 
of heterochrony identified by our approach may be asso- 
ciated with specific biological functions. We first tested 
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the characteristics of genes in the three time shift clus- 
ters using published data for gene expression specific to 
human brain gray and white matter [21], as well as on 
cell type-specific expression in the mouse central ner- 
vous system [22]. Genes in Cluster 1 were significantly 
enriched among genes preferentially expressed in 
human gray matter (hypergeometric test [HT] p = 
0.001), while genes in the other two clusters did not 
show any cell type specificity ([HT] p > 0.5) (Figure 4E). 
Consistently, genes in Cluster 1, but not in the other two 
clusters, tended to be enriched among neuron-specific 
genes ([HT] p = 0.058) (Figure 4E). These results indi- 
cate that differences in time shift patterns, identified by 
DTW-S, might reflect differences in ontogenetic timing 
of gene expression changes among different human 
brain cell types, and between histological locations. 
Development of the human and macaque prefrontal cortex 
and cerebellum 

To further test the utility of the DTW-S in uncovering 
novel biological features, we studied time shift patterns 
in human and rhesus macaque brain development in 
two distinct regions: prefrontal cortex and cerebellum. 
For this, we used a published expression dataset con- 
taining time series measured in the dorsolateral prefron- 
tal cortexof 23 humans and 26 rhesus macaques, and in 
the cerebellum of 22 humans and 24 rhesus macaques 
[23]. In both species, the dataset covered most of the 
lifespan, humans: 0-98 years, macaques: 0-28 years 
(Additional file 1, Figure SIB). 

Based on these data, we selected 1084 genes in the 
cortex and 950 genes in the cerebellum that satisfied the 
following criteria: (1) significant expression changes with 
age (using polynomial regression models, at F-test p < 
0.05), (2) significant expression differences between 
humans and rhesus macaques (using analysis of covar- 
iance (ANCOVA), at F-test p < 0.05), (3) significant 
positive correlation between human and rhesus macaque 
expression profiles (Pearson correlation p < 0.05 and r > 
0) and, (4) significant time shift between human and 
rhesus macaque expression profiles (FPR = 11.7% for 
cortex and 11.5% for cerebellum). Taking the union of 
the genes that passed these criteria in the cortex or cer- 
ebellum yielded 1735 genes. 

For these genes, we determined the direction of the 
time shift (acceleration, neoteny, or none) between 
human and rhesus macaque time series, in both cortex 
and cerebellum, by aligning rhesus macaque data to the 
human age scale. If at least 70% of time points showed 
consistent time shift direction, we classified this gene as 
accelerated or neotenic. Following this procedure, the 
vast majority of genes were classified as neotenic, with 
1444 (83%) and 1469 (85%) in the cortex and cerebel- 
lum respectively. This result is consistent with the faster 



rate of macaque brain development and maturation that 
has been previously reported [24]. 

Overall, time shift measurements correlated positively 
between the prefrontal cortex and cerebellum (Additional 
file 1, Figure S8), with 560 out of 1750 genes tested 
showed strong positive correlation (Pearson p < 0.05, r > 
0.5). The time shift profiles of these genes could be 
assigned to four clusters using the k-means clustering 
method (Figure 5A). Notably, the vast majority of these 
genes (452 genes in clusters 1 and 2) showed nearly iden- 
tical time shift profiles in the prefrontal cortex and cere- 
bellum. Does this similarity of time shift profiles reflect 
similarity in gene expression profiles in the two brain 
regions? Indeed, for some genes in both clusters 1 and 2, 
gene expression profiles in cortex and cerebellum fol- 
lowed the same trajectories (Figure 5B). Other genes, 
however, showed clearly distinct expression profiles in 
the two brain regions, while still sharing the same time 
shift profile. In a search for a biological explanation of 
this phenomenon, we tested cell type specificity and 
enrichment in biological functions specified by Gene 
Ontology (GO) annotation [25]. We found that genes 
grouped into Clusterl based on their time shift profiles, 
and showing identical expression patterns in cortex and 
cerebellum (CU-groupl), were enriched in white matter 
([HT] p = 0.0002), and were annotated in GO terms: 
"lipid metabolic process" ([HT] p = 0.007) and "cellular 
lipid metabolic process" ([HT] p = 0.007). By contrast, 
Clusterl genes showing different expression patterns in 
cortex and cerebellum (Cll-group2) were enriched in 
gray matter ([HT] p = 0.0004) and annotated by GO as 
"nervous system development" ([HT] p = 0.05) (Figure 
5C). Similarly, genes grouped into Cluster2 based on 
their time shift profiles, and showing identical expression 
patterns in cortex and cerebellum (Cl2-group2), were 
enriched in white matter ([HT] p = 0.02) and mature oli- 
godendrocytes ([HT] p = 0.03), while genes showing dif- 
ferent expression profiles (Cl2-groupl) were enriched in 
gray matter ([HT] p = 0.006). Thus, our results revealed 
an interesting biological phenomenon: within one species, 
ontogenetic profiles are shared between the prefrontal 
cortex and cerebellum for genes expressed in white mat- 
ter, but distinct for genes expressed in gray matter. 
Importantly, however, the time shift between the human 
and rhesus macaque ontogenetic profiles is perfectly syn- 
chronized for both white and gray matter genes. On an 
organismal level this observation might not be surprising. 
Changes in the rate of ontogenesis might be expected to 
operate on the brain as a whole, leading to synchronized 
delay in white and gray matter development in humans, 
compared to rhesus macaques. Our results confirm that, 
on the gene expression level, such synchronization can 
indeed be observed. 
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Figure 5 Human-rhesus macaque time shift analyzed by DTW-S. (A) Four patterns of time shift between human and rhesus macaque 
expression time series, showing high time shift profile correlation (r: 0.5-1) between cortex and cerebellum. The time shift was calculated by 
aligning rhesus expression time series to human time series. The numbers above the graphs represent numbers of genes within each cluster. 
The graphs show mean time shift of all genes within a cluster (points) and the time shift interval for individual genes within a cluster (5%-95%). 
Positive values of time shift indicate that human ages are mapped to younger ages in the rhesus macaque time series, indicating slower/delayed 
human development. Time shift shown on the y-axis is calculated as the base-two logarithm ratio between rhesus macaque age and human 
age at aligned expression time points. (B) Patterns of expression profiles of human genes in the time shift Cluster! and Cluster2 shown in (A). 
The numbers above each graph represent genes within a group. The graphs show mean expression of all genes within a group (points) and 
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expression level variance (error bars) for individual genes within a cluster (5%-95%). (C) Properties of genes in the time shift Cluster! and 
Cluster2. The left and the right panels show the proportions of genes in each of the four groups of Cluster!, and Cluster2 respectively, that 
overlap with genes with expression specific to gray matter, white matter or non-specific. 



Conclusions 

Comparisons of developmental patterns across closely 
related species are playing an increasingly important role 
in extracting meaningful information from biological time 
series. By modifying the dynamic time warping algorithm 
(DTW), we have designed an effective tool for time series 
alignment (DTW-S). Our simulation results show that this 
method is effective in calculating time shift between two 
time series, even when the proportion of noise is 20-30% 
of the total variance. Furthermore, this method performs 
well for expression profiles containing both recurrent and 
non-recurrent changes, and can estimate variation in the 
amplitude and direction of the time shift. 

When we applied our method to a published gene 
expression dataset of human, chimpanzee and rhesus 
macaque brain development and maturation, we 
obtained robust and reproducible time shift estimates 
consistent with previous observations [15]. Furthermore, 
our method allowed us to classify genes into distinct 
categories according to their time shift patterns. This 
provided additional insights into biological mechanisms 
underlying human-specific brain development and 
maturation, which could not be deduced from the gene 
expression profile data alone. 

Applying our method to a gene expression dataset of 
human and rhesus macaque brain development and aging, 
we found that genes showing synchronized time shift 
between the species, in the prefrontal cortex and cerebel- 
lum, do not always follow the same expression profiles in 
the two brain regions. Notably, genes showing both syn- 
chronized time shift between human and macaque onto- 
genetic trajectories, and synchronized expression patterns 
in the prefrontal cortex and cerebellum, were preferen- 
tially expressed in brain white matter. By contrast, genes 
showing synchronized time shift between human and 
macaque ontogenetic trajectories, but different expression 
patterns in the prefrontal cortex and cerebellum, were pre- 
ferentially expressed in brain gray matter. 

Taken together, these two examples demonstrate that 
the combination of gene expression time series profiles 
with ontogenetic time shift estimates provides additional 
information revealing the biological properties of the 
investigated system. The development of DTW-S algo- 
rithm, freely available as the R package "TimeShift", 
should facilitate the application of this approach to 
further studies. 



Additional material 



Additional file 1: Supplementary figures, includes all additional figures 
mentioned in this paper. 
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